library(MASS)

simulations <- 1000
sampleSizes <- c(250,500)
treatmentDifferences <- c(0.05,0.02,0.01)
powerAnalysis <- data.frame(n=integer(),
                            treatmentDiff = integer(),
                            pValue = integer()
                            )

#power analysis for support of draconian measures
for(n in sampleSizes)
{
  for(treatmentDiff in treatmentDifferences)
  {
    for(i in 1:simulations)
    {
      #r <- sqrt(treatmentDiff)
      X1 <- rnorm(n, mean = 0, sd = 1)
      #X2 <- rnorm(n, mean = 0, sd = 1)
      #X3 <- r*X1+sqrt(1-r^2)*X2
      treatmentA <- X1
      #treatmentB <- X3
      #the values of the mean have been found through simulation
      treatmentB <- rnorm(n, mean = ifelse(treatmentDiff == 0.05, 0.46, ifelse(treatmentDiff == 0.02,0.286,0.2)), sd = 1)
      simData <- data.frame( DV = c(treatmentA,treatmentB),
                             IV = c(rep("A",n),rep("B",n)))
      s <- summary(lm(DV~IV, data = simData))
      powerAnalysis[nrow(powerAnalysis)+ 1,] <- list(n, treatmentDiff, s$coefficients[2,4])
    }
  }
}
powerAnalysis$truePos <- powerAnalysis$pValue <= 0.05
aggregate(powerAnalysis$truePos, 
          by = list(powerAnalysis$treatmentDiff,powerAnalysis$n),
          FUN = mean)


#power analysis for donation decision.
#We expect that 50% of subjects will donate
powerAnalysis <- data.frame(n=integer(),
                            treatmentDiff = integer(),
                            pValue = integer()
                            )
for(n in sampleSizes)
{
  for(treatmentDiff in treatmentDifferences)
  {
    for(i in 1:simulations)
    {
      treatmentA <- rbinom(n, size = 1, prob = 0.5)
      treatmentB <- rbinom(n, size = 1, prob = ifelse(treatmentDiff == 0.05, 0.28, ifelse(treatmentDiff == 0.02, 0.36,0.4)))
      simData <- data.frame( DV = c(treatmentA,treatmentB),
                             IV = c(rep("A",n),rep("B",n)))
      m <- glm(DV~IV, data = simData,family=binomial(link='logit'))
      s <- summary(m)
      s
      #rsq(m)
      powerAnalysis[nrow(powerAnalysis)+ 1,] <- list(n, treatmentDiff, s$coefficients[2,4])
    }
  }
}
powerAnalysis$truePos <- powerAnalysis$pValue <= 0.05
aggregate(powerAnalysis$truePos, 
          by = list(powerAnalysis$treatmentDiff,powerAnalysis$n),
          FUN = mean)




